Reservoir simulation is widely used in the petroleum industry to analyze and forecast the behavior of fluid flow in hydrocarbon bearing subterranean reservoirs, such as gas or oil fields. For instance, simulation inputs can be varied for a reservoir model to project how certain changes might influence a reservoir's future performance. The simulation results are typically used for making reservoir management decisions, such as whether to develop new reservoir fields or how to efficiently manage production of current fields. Developing accurate and more computationally efficient methods for use in reservoir simulators has become increasingly important in the petroleum industry.
In general, reservoir simulation involves the numerical solution of a system of equations that describes the physics governing the complex behaviors of multi-component, multiphase fluid flow in the naturally porous media of a subterranean reservoir. The system of the equations is typically in the form of coupled nonlinear partial differential equations (PDE's). The PDE's are discretized in time and space on a given grid, and discrete equations are solved, such as by an iterative process, for a series of time steps until a prescribed time is reached. At each time step, linearization of the nonlinear system of equations (e.g., Jacobian construction), solving the linear system, and computing a subsequent system of equations is performed. A display representing the fluid flow in the reservoir being modeled can be generated using the simulation results at each time step.
For example, a grid (structured or unstructured) can be imposed upon an area of interest in a reservoir model to define a plurality of cells, each having one or more unknown properties associated therewith. Examples of unknown properties can include, but are not limited to, fluid properties such as pressure, temperature or water saturation, and rock properties such as permeability or porosity. A matrix can be constructed to represent the gridded region of interest and to solve for the unknown variables. For a grid of m cells and ki unknowns per cell i, i=1, . . . m, giving a total of n unknowns, the linear systemAx=b  Equation (1)is obtained by discretization of the partial differential equations describing flow within the reservoir and wells. Here, A is a known square n×n matrix, x is an n-dimensional vector of the variables representing unknown properties of the cells to be found by solving the system of equations, and b is an n-dimensional vector of known quantities. Block vector x and block vector b are the same length. If a fully implicit time discretization is used, the matrix will generally be constant block-size. If adaptive implicit (AIM) time-stepping is used, such as IMPSAT (implicit pressure and saturation) or IMPES (implicit pressure explicit saturation), the matrix can be variable block-size.
The unknowns are typically ordered cell by cell. Matrix A can be partitioned by cells yielding the blocked matrixA=[aij]i,j=1, . . . m  Equation (2)The nonzero aij sub-matrices correspond to the cell neighbors of cell i determined by the finite difference stencil and the block aij is ki×kj for i, j=1, . . . , m.
The solution of the linear system of equations n be a very computationally-intensive task. Generally, linear solvers utilize either direct methods or iterative methods to determine the solution. In the direct method, Gaussian elimination is used such that matrix A is factorized into a product of lower triangular matrix, L, and upper triangular matrix, U. Accordingly, A=LU in the direct method. However, for large sparse matrices, computation of triangular matrices L and U becomes prohibitively expensive as the number of non-zero elements in each factor is very large. In iterative methods, the linear system is solved using approximations to matrix A. For example, an incomplete lower-upper (ILU) factorization can be used, instead of a full factorization as in the direct method. Here, a product of sparse factors L and U are computed such that their product approximates matrix A (A≈LU). When employing an iterative method, the solution is updated repetitively until convergence is reached. Unfortunately, standard iterative methods converge very slowly for large systems of linear equations because the number of iterations typically increases as the number of unknowns increases.
Preconditioning can be used to decrease the number of iterations used by the iterative method when solving the linear system. In preconditioning, matrix A of Equation (1) is multiplied by a preconditioning matrix, often called a “preconditioner” for brevity, such that the linear system is more suitable for numerical solution. For example, matrix A can be transformed into a right preconditioned system, given as AM−1Mx=b, or a left preconditioned system, given as M1(Ax−b)=0. Here M1 represents an inverse of preconditioning matrix M. In either case, the preconditioner can be used to increase the rate of convergence when solving Equation (1) since the number of iterations typically grows more slowly, or not at all, with the size of the problem (i.e., the number of unknowns).
Single or multi-stage preconditioners (AM−1 or M−1A) are frequently combined with Krylov subspace methods for solving the fully implicit or adaptive implicit systems of equations. Examples of well-known iterative procedures used to accelerate convergence in a Krylov subspace include GMRES, which is taught by Y. Saad and M. H. Schultz in “GMRES: A Generalized Minimal Residual Algorithm for Solving Non-symmetric Linear Systems,” SIAM J Sci. Stat. Comp., 7(3): 856-869, 1986. One skilled in the art will appreciate that other Krylov subspace methods and algorithms can also be used to accelerate domain decompositions such as the Flexible Generalized Minimal Residual Algorithm (FGMRES) taught by Y. Saad in “A Flexible Inner-Outer Preconditioned GMRES Algorithm,” SIAM J. Sci. Stat. Comp., 14(2): 461-469, 1993; the Orthogonal Minimization (ORTHOMIN) algorithm as taught by P. K. W. Vinsome in “Orthomin, an Iterative Method for Solving Sparse Sets of Simultaneous Linear Equations,” which was presented at the SPE Symposium on Numerical Simulation of Reservoir Performance, Feb. 19-20, 1976 in Los Angeles, Calif.; the Bi-Conjugate Gradient Stabilized (Bi-CGSTAB) algorithm as taught by H. A. van der Vorst in “Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems,” SIAM J. Sci. Stat. Comp., 13(2): 631-644, 1992; and the Generalized Conjugate Residual (GCR) algorithm as taught by S. C. Eisenstat, H. C. Elman, and M. H. Schultz in “Variational Iterative Methods for Nonsymmetric Systems of Linear Equations,” SIAM Journal of Numerical Analysis, 20(2) 345-357, 1983.
A commonly used preconditioner in reservoir simulation is the Constrained Pressure Residual (CPR) preconditioner, which is taught by John R. Wallis in “Incomplete Gaussian Elimination for Preconditioning of Generalized Conjugate Gradient Acceleration,” which was presented at the SPE Reservoir Simulation Symposium, Nov. 15-18, 1983 in San Francisco, Calif. (SPE 12265). The CPR preconditioner exploits the characteristics of pressure being elliptic (or substantially elliptic) and a “stiff” variable by constructing a two-stage preconditioner comprising a first stage pressure matrix App=WTAC and second stage implicit preconditioning matrix M.MCPR−1=M−1[I−ACApp−1WT]+CApp−1WT  Equation (3)Here, WT is an m×n restrictor matrix and C is an n×m in interpolation matrix. Examples of methods used to determine these matrices include implicit pressure explicit saturation (IMPES) formulations, such as True-IMPES and Quasi-IMPES. Additional details on IMPES reduction methods can be found in U.S. Pat. No. 7,516,056. Such methods approximately decouple the pressure variables from the non-pressure variables. The pressure solutions in CPR are usually accomplished by one or more steps of preconditioned GMRES, FGMRES, ORTHOMIN, BICGSTAB, or simply one preconditioning step. Examples of pressure preconditioners include Algebraic Multi-Grid (AMG) methods, and incomplete lower-upper (ILU) triangular factorizations including ILU(k), ILUT(τ), and ILUT(p,τ) factorizations. The value k denotes the level of infill allowed for ILU(k). The value τ denotes the infill drop tolerance. Similarly in ILUT(p,τ) only the p largest magnitude elements in the L part of each row are kept, which are also greater than the drop tolerance. The same criterion is applied to each row of U and the diagonal element is always kept. Accordingly, ILUT (incomplete lower-upper with threshold) factorization computes the entire row of the ˜L and ˜U matrices, and maintains values greater than a predetermined threshold such that the amount of fill-in is controlled during factorization by dropping entries smaller than a prescribed threshold.
The CPR method is very effective in improving convergence over standard single-stage ILU(k) and ILUT methods. The second stage preconditioning of the CPR method is typically ILU(k) on the multi-variable per cell block matrix where the value of k is zero (0) for easy problems and one (1) for more difficult problems. However, the implementation of block ILU(k) does not exploit sparsity within the sub-matrix blocks and does not use the magnitude of the infill terms to adaptively modify the factorization sparsity pattern to improve accuracy. The ILUT methods are typically used for unblocked matrices, and are useful in addressing accuracy by its adaptive sparsity pattern. However, they typically cannot be applied directly to highly non-symmetric fully implicit or adaptive implicit matrices without using pivoting, which is typically computationally expensive.
While the above methods can be used in linear solvers, a robust and computationally efficient method is needed for solving linear systems of equations that avoids the aforementioned shortcomings, particularly when being applied to large scale reservoir simulation problems. Moreover, any new method should also produce excellent convergence rates.